package twoD.hofem;

import inf.jlinalg.*;

public class FMath {
	// Adds two functions
	public static FunctionRToR add(final FunctionRToR u, final FunctionRToR v) {

		checkCompatibility(u, v);

		return new FunctionRToR() {
			public double derivativeAt(double x) {
				return u.derivativeAt(x) + v.derivativeAt(x);

			}

			// Breakpoints of the resultant function is union of the breakpoints
			// of the input functions
			public double[] getBreakpoints() {
				double[] a = u.getBreakpoints();
				double[] b = v.getBreakpoints();
				return combine(a, b);

			}

			public Interval getDomain() {
				return u.getDomain();
			}

			public int getP() {
				int a = -1;
				if (u.getP() != -1 && v.getP() != -1) {
					a = Math.max(u.getP(), v.getP());
				}
				return a;
			}

			public double valueAt(double x) {
				return u.valueAt(x) + v.valueAt(x);
			}
		};

	}

	// Checks if the two functions have the same domain
	private static void checkCompatibility(final FunctionRToR u,
			final FunctionRToR v) {
		Interval du = u.getDomain();
		Interval dv = v.getDomain();
		double lb = Math.abs(du.getLower() - dv.getLower());
		double ub = Math.abs(du.getUpper() - dv.getUpper());
		if (Math.max(lb, ub) > 1e-12) {
			throw new IllegalArgumentException("Domains incompatible");
		}
	}

	// Multiplies two functions
	public static FunctionRToR multiply(final FunctionRToR u,
			final FunctionRToR v) {
		checkCompatibility(u, v);

		return new FunctionRToR() {
			public double derivativeAt(double x) {
				return u.valueAt(x) * v.derivativeAt(x) + u.derivativeAt(x)
						* v.valueAt(x);
			}

			// Breakpoints of the resultant function is union of the breakpoints
			// of the input functions
			public double[] getBreakpoints() {
				double[] a = u.getBreakpoints();
				double[] b = v.getBreakpoints();
				return combine(a, b);

			}

			public Interval getDomain() {
				return u.getDomain();
			}

			public int getP() {
				if (u.getP() == -1 || v.getP() == -1) {
					// in some cases it is not true, e.g, u= e^x and v= x/e^x ,
					// u*v= x is a polynomial
					return -1;
				} else {
					return u.getP() + v.getP();
				}

			}

			public double valueAt(double x) {
				return u.valueAt(x) * v.valueAt(x);
			}
		};
	}

	// Multiplies a function with a scalar value
	public static FunctionRToR multiply(final double a, final FunctionRToR v) {

		return new FunctionRToR() {
			public double derivativeAt(double x) {
				return a * v.derivativeAt(x);
			}

			// Breakpoints of the resultant function is the same as the input
			// function
			public double[] getBreakpoints() {
				return v.getBreakpoints();
			}

			public Interval getDomain() {
				return v.getDomain();
			}

			public int getP() {
				return v.getP();
			}

			public double valueAt(double x) {
				return a * v.valueAt(x);
			}
		};
	}

	// Composes two functions : (g o f) ( x ) = g( f ( x ) )
	public static FunctionRToR compose(final FunctionRToR u,
			final FunctionRToR v) {

		return new FunctionRToR() {
			public double derivativeAt(double x) {
				return u.derivativeAt(v.valueAt(x)) * v.derivativeAt(x);

			}

			public double[] getBreakpoints() {
				return u.getBreakpoints();

			}

			public Interval getDomain() {
				return v.getDomain();
			}

			public int getP() {
				if (u.getP() == -1 || v.getP() == -1) {

					return -1;
				} else {
					return u.getP() * v.getP();
				}

			}

			public double valueAt(double x) {
				return u.valueAt(v.valueAt(x));
			}
		};

	}

	// In order to sort the elements of an array

	public static double[] Sort(double[] a) {
		double temp;
		int pos;
		for (int i = 1; i < a.length; i++) {

			temp = a[i];
			pos = i - 1;

			while ((pos >= 0) && (temp < a[pos])) {
				a[pos + 1] = a[pos];
				pos--;
			}
			a[pos + 1] = temp;
		}
		return a;
	}

	// In order to sort the elements of an array,which combines two arrays,
	// in ascending order by avoiding repetition of the same numbers in the new
	// array

	public static double[] combine(double[] a, double[] b) {
		if (a == null) {
			return b;
		} else if (b == null) {
			return a;
		} else {
			int k = a.length + b.length;
			for (int i = 0; i < a.length; i++) {
				for (int j = 0; j < b.length; j++) {
					if (a[i] == b[j]) {
						k--;
					}
				}
			}
			double[] d = new double[a.length + b.length];
			for (int p = 0; p < a.length; p++) {
				d[p] = a[p];
			}
			for (int q = 0; q < b.length; q++) {
				d[q + a.length] = b[q];
			}
			d = Sort(d);
			double[] c = new double[k];
			int i = 0;
			int f = 0;
			while (i < d.length - 1) {
				if (d[i] != d[i + 1]) {
					c[f] = d[i];
					if (i == d.length - 2) {
						c[f + 1] = d[i + 1];
					}
					i = i + 1;
				} else {
					c[f] = d[i];

					i = i + 2;
				}
				f++;
			}
			return c;
		}

	}

	/**
	 * This method determines the coefficients c[] to fit the given function U
	 * from the given basis functions b Uh = c1*b1 + c2*b2...
	 * 
	 * @param u
	 *            FunctionRToR
	 * @param b
	 *            FunctionBasisOnR
	 * @return Coefficients matrix c[]
	 * 
	 *         Author: Team4
	 */
	public static double[] fit(FunctionRToR u, FunctionBasisOnR b) {
		FunctionRToR[] d = b.getBasisFunctions();
		int n = d.length;
		BilinearFormOnRToR bf = new L2ScalarProduct();
		double[][] A = new double[n][n];
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = bf.valueAt(d[i], d[j]);
			}
		}
		double[] f = new double[n];
		for (int i = 0; i < n; i++) {
			f[i] = bf.valueAt(u, d[i]);
		}
		// Solving system linear equation A*c = f
		ILSESolver solver = new GeneralMatrixLSESolver();
		QuadraticMatrixInfo aInfo = solver.getAInfo();
		IMatrix a = solver.getA();
		aInfo.setSize(n);
		solver.initialize();
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				a.set(i, j, A[i][j]);
			}
		}
		try {
			solver.solve(f);
		} catch (SolveFailedException e) {
			System.out.println("Solve failed : " + e.getMessage());
		}
		return f;

	}

	public static FunctionR2ToR multiply(final FunctionR2ToR u,
			final FunctionR2ToR v) {
		return new FunctionR2ToR() {

			public double[] gradientAt(double xi1, double xi2) {
				throw new Error();
			}

			public double valueAt(double xi1, double xi2) {
				return u.valueAt(xi1, xi2) * v.valueAt(xi1, xi2);
			}
		};
	}
}
